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^N ■ The method used earher for analysis of correlated nanoscopic systems is extended to infinite 

(periodic) s-band like systems described by the Hubbard model and its extensions. The optimized 
single-particle wave functions contained in the parameters of the Hubbard model (the hopping t and 
the magnitude of the intraatomic interaction U) are determined explicitly in the correlated state for 

VO ' the electronic systems of various symmetries and dimensions: Hubbard chain, square and triangular 

planar lattices, and the three cubic lattices (SC, BCC, FCC). In effect, the evolution of the electronic 
properties as a function of interatomic distance R is obtained. The model parameters in most cases 
do not scale linearly with the lattice spacing and hence, their solution as a function of microscopic 
D I parameters reflects only qualitatively the system evolution. Also, the atomic energy changes with 

< I . R and therefore should be included in the model analysis. The solutions in one dimension (D = 

1) can be analyzed both rigorously (by making use of the Lieb-Wu solution) and compared with 
the approximate Gutzwiller treatment. In higher dimensions {D — 2, 3) only the latter approach 

T^ ■ is possible to implement within the scheme. The renormalized single particle wave functions are 

a almost independent of the choice of the scheme selected to diagonalize the Hamiltonian in the Fock 

space in D = 1 case. The method can be extended to other approximation schemes as stressed at 
1-^ , the end. 

ti . PACS numbers: 71.27.+a 71.30.+h, Tl.lOFd 

I. INTRODUCTION 

> 

^\\ ' The question of combining in an explicit manner inter-electronic correlations with single-particle (band) calculations 

■^ , is very important for the systems for which Coulomb interaction between electrons is comparable to or even larger than 
'nJ" ' the kinetic (bare band) energy of electronaii^. With respect to this, methods, starting from band calculations, such 
^N ^ as LDA-t-U"^ or LDA-f DMFT^ have been devised and they work well for quite few systems. The methods provide e.g. 
[ — ■ the photoemission spectrum^, the overall band splitting at the Fermi surface (i.e. the stability of the Mott insulating 
^^ ' state) optical spectrum, etc. The band structure calculations allow also for estimation of the Hubbard-interaction 
OO . parameter U, when the Wannnier functions are determined first^. In all these methods the question of counting twice 
the electron-electron Coulomb extended interaction arises when the effective single-particle potential contains them. 
A systematic approach base on taking into account Coulomb interactions between electrons in an exact manner 
first and only then determining the renormalized single-particle wave functions contained in the model parameters by 
. , a proper Euler variational procedure. In that situation we allow for an adjustment of single-particle wave functions 
j^ ' in the correlated-electron state and only then calculate them explicitly. Such a reverse method called EDABI (Exact 
Diagonalization with Ab Inition Aproach) has been devised and employed to nanoscopic chains and clusters^. It 
provides the evolution of the correlated-system properties as a function of interatomic distance R, not only as a 
function of model parameters, which are difficult to be measure. For example^, EDABI provides new results such 
as e.g. a systematic evolution of the statistical distribution function as a function of increasing R (from Fermi-like 
function to a continuous momentum distribution reflecting electron localization on parent atoms) or a magnetic Slater- 
like splitting of the electronic states without the appearance of long-range antifcrromagnctic ordering^. Needless to 
say that this method avoids in an explicit manner counting twice the interaction between the particles. 

The purpose of this paper is to generalize and test the EDABI-type approach for extended systems of arbitrary 
dimension {D = 1,2,3) described by the (parametrized) extended Hubbard model. Only in the D = 1 case it is 
possible to compare an exact (Lieb-Wu, LW) solution with the approximate Gutzwiller- wave- functions (GWF), and 
the Gutzwiller-ansatz (GA) solutions. In higher dimensions, we calculate the single-particle properties starting from 
GA. What is surprising, at least for Z? = 1, is the relative insensitivity of the detailed shape of the renormalized (by 
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FIG. 1: Flowchart of the approach combining exact diagonalization in the Fock space followed by the single-particle wave- 
function optimization in the Hilbert space for a correlated state. For details see main text. 

correlations) single-particle wave function to the method selected to diagonalize the many-particle Hamiltonian in the 
Fock space. In general, our method of approach completes the solution of the parametrized models in the sense that 
it yields the evolution of the correlated many-particle systems properties as a function of interatomic distance, as well 
as provides the shape of the Wannier functions in the correlated state. 

The structure of the paper is as follows. In Section II and III we overview briefly our method, whereas in Section 
IV the extended Hubbard-chain properties are analysed in detail. In Sections V and VI selected two- and three- 
dimensional lattices are considered respectively. Section VII contains a brief discussion and an overview. In Appendices 
A and B we provide some formal details of the calculations. 

II. PARAMETERIZED MODELS SUPPLEMENTED WITH THE SINGLE-PARTICLE BASIS 
OPTIMIZATION: A BRIEF SUMMARY OF EDABI 



Before placing our work in the literature of the subject, we first characterize the essence of our approach^ from 
the formal side. The method of the optimized single-particle wave functions incorporates the first and the second 
quantization schemes. Namely, if ^(ri, ..., rjv) describes the A'^-particle wave function in the Schrodinger representation 
(i.e. in Hilbert space), then in the second-quantization representation (i.e. in the Fock space) this state can be 
represented by^ 



1$) --L=yd3ri...dV^(ri,...,rjv)^t(i.^)...^t(r^)|o), 



(1) 



where, ^^(rj) = f 'i't(ri)'4'!(ri) j, denotes the field operators representing the particle creation at point r, with 

spin a =1,1= ±1. Utilizing the anticommutation relations for the field operators one can easily find the inverse 
representation which has the form 



*(ri,...,rAr) = 



1 



/M 



^(ri)...^(rjv)|$) 



(2) 



So, the two schemes are equivalent if only the anticommutation relations of the field operators are defined. Here, 
the field operator is defined in the single-particle basis of Wannier functions Wi (r) as follows 



i=l,o-=±l 
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(3) 



where the summation over i runs over a complete basis of single-particle states. The basis {wi(r)} can be arbitrary, 
provided it is complete in the quantum-mechanical sense. In most situations we introduce a model, i.e. select the 
subset of a complete basis {{]. For example we replace this complete basis by a finite subset {u)i(r)}, i — 1, ..., M of 
the functions, connected directly to the problem at hand. Hence, in the single narrow-band situation we have thatji 



M 



^(r)«^z^.(r)a„ (4) 



where in what follows we set M = N and {wi(r)} is the single-particle- Wannier- function basis (to be defined later). 
In effect, the approximate states in the Fock space are defined through 



M 



^ E Q....,«4-<Jo)' (5) 



and the iV-particle wave function in this approximate basis has now the form 

M M 

*(ri,...,rjv) = ^ E E (0|aiiv-anajt--aj>|0> xC'ji-iivW'.i(ri)...Wi„(rAr). (6) 

The coefficients Cj-^,,,jj^ can be calculated by either direct Hamiltonian diagonalization or by Lanczos method for fi- 
nite systems'^'^, whereas the normalized (optimized) wave functions {wi{r)} are determined by the procedure described 
below. Note that the second-quantization formalism separates the many-particle function aspect of the problem, which 
is contained in the coefficients C'ji...j„, from the wave-mechanics aspect of determining the basis {wi(r)}. In fact, the 
two are intertwinned. Namely, we perform the diagonalization of the second-quantized Hamiltonian for selected (and 
fixed) single-particle basis first and only then optimize the basis {wi(r)} with the help of a variational approach. 

We now describe how to combine the second and the first quantization schemes on the example of extended Hubbard 
model. The general Hamiltonian for interacting particles is 

-H" = E*y"tajT + :5 E Vi3kiata+^,aia,aka, (7) 
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where 



U, = {w^Hi \w,) ~ / d'vw*{v)H^Wj(v), (8) 



is the hopping integral with 



being the Hamiltonian of single bare electron (a.u.=atomic units), and is the Hamiltonian for single electron in the 
system under consideration, and 



Vijkl = {w,Wj\V \wkWl) = / d-^Vid'^Y2W*{vi)w*{r2)V{Yi - r2)wk{Yl)wi{Y2) (10) 

is the amplitude of classical Coulomb interaction rewritten in the first-quantization language. Note that here i labels 
complete set of quantum numbers except spin. The ground state energy is then determined by the expression 

Eg = {H)=Y^ Uj U\^aj„\ + 7^ Yl Viju (atapaia'aka), (H) 



2 

ij(T ijklaa' 



where the ground-state averages (...) are determined for the many particle ground state |$) = |$g) for fixed {wi{r)}. 
This means, that the trial basis {^^(r)} i = l,...,Af entering tij and Vijki must be optimized, i.e. the ground 
state energy must be a minimum within the class of trial (incomplete) basis of wave-functions. Obviously, such an 
optimization would not be necessary if the basis were complete. In the present approach this means that we have 
to construct the functional variational scheme to determine the optimized trial basis {wi{r)}. This formal procedure 
provides us with the evolution of the system properties as a function of the lattice parameter R. The method is 
summarized schematically in Fig. [TJ 

III. THE MODEL 
A. Extended Hubbard model 

To describe a single-band model of interacting fermions we start with the extended Hubbard Hamiltonian 

H = ea'^ni+ ^ tijata^^jcr + U^ ni-^riii + ^ KijUiUj + ^ V^on(Rj - ^i), (12) 

where now i labels a Wannier orbital centered at i = R^, Eq = ta^ is the atomic energy per site, tij - the hopping 
integral between the sites i and j ^ i, t/ is the magnitude of intrasite Coulomb interaction, whereas Kij is the 
corresponding quantity for electrons located on sites i and j with j ^ i. Finally, 

Vion-ion — p5 5— I = -f^ (13) 

|iv,i — -Ttjl -Ky- 

is the classical Coulomb interactions between the cations located at the sites i and j in atomic units (a.u.). One 
should note that above i and j mean the ionic positions, so the Hubbard model already does not base on a complete 
set of the single-particle wave-functions {wi(r)}, i.e. neglects p, d, /, etc. Wannier states. 

Even though we study here the Hubbard model, we have to include the intersite Coulomb interaction if the atomic 
limit is to be recovered properly in the limit of interatomic distance Rtj — > 00. In order to achieve that, we first 
represent the intersite term in the form — 



i<j i<j i<j 

where N^ = J2 ''^i i^ ^^^ total number of electrons, N is the number of lattice sites, and Srii = n^ — 1. In the Mott 

i 

insulating state the number of electrons N^ = N i.e. in the half-filled case, {Sui) = 0. This condition defines the 
Mott-Hubbard state. In this state, we have that 

y^ KijUiUj = ^ Kij. (15) 

i<j i<j 

In effect, Hamiltonian IT^ in that limit assumes the form 

^ = Ea^E"* + E ^*i««T«J'^ + [/^n,|n,|, (16) 



i<j,(y 



where 



Cfj On 



'^4e(^. + i|-)> (17) 



is the effective atomic energy which reduces to the true atomic energy in the hmit of large interatomic separation. 
The inclusion of this energy part assure the reduction of e^J^ to that of the isolated atoms for i?y ^ oo. In that 
situation, e^J^ is not a constant quantity, as would be the case for a parameterized model. The ground state energy 
is 

One sees that in this expression we have neglected the contribution of charge fluctuations ^ "^ Kij {SuiSrij). 
Obviously, it vanishes in the Hartree-Fock approxiamtion. It constitutes the only approximation, under which the 
Lieb-Wu solution can be applied to the extended Hubbard chain as well. 

B. Constructions of the single-particle basis and the construction of the attractive periodic potential 

We first construct the single-particle wave-function basis {wi{r)} entering the expressions for t, U, and Kij in the 
tight binding approximation i.e. assume that 

z 

u;,(r)=/3vl',(r)-7^vl/,(r), (19) 

where z is the number of nearest neighbors, j3 and 7 are the mixing coefficients, and 



^^(r) = W—e-"!---^-! (20) 

V TT 

is the Is atomic wave function centered on the site i. In concrete calculations, they are represented by the Gaussians, 
as defined below. The parameters (3 and 7 are selected to fulfill the atomic-limit properties 



where ag is the atomic Bohr radius and 



lim a = ctQ = l/flOi (21a) 

-R— ►00 



lim /3 = 1, lim 7 = 0. (21b) 

R — >oo R — *oo 



As the functions {wi(r)} are orthogonalized atomic orbitals, we have that 

< Wi|wi >= 1, <Wi\wj>^0. (22) 

This conditions lead to the explicit expressions of the form 



. A + ^W^BSl 

[2A2 - BSi - zASl + 2{A - zSly A^ - B5i]i/2 ' ^ > 



and 

Si 



[2A2 - BSi - zASl + 2{A - zSf)y'A^ - BSi]^/^ ' 
with 



(24) 



and 



A^ J2 (*^l%w)= E (*«l*.ia))' (25) 
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FIG. 2: Effective "periodic" attractive Coulomb potential seen by electron localized on site j of a linear chain of (+e) kations 
and composed of: (a) 6 wells (j e 52(i) U S2{i + 1)) and (b) 22 wells {j G 5*10(1) U S\o{i + 1)). For details see main text. 



whereas 



is the overlap integral. In general, 



5i = (*rJ*r,^,) 



5„ = (*rJ*r,, 



(27) 



(28) 



is the overlap that with the n-th coordination sphere. The symbol j(i) in the above equations labels j neighbors of 
the site i. The explicit form of the coefficients A and B for the structures studied is provided in Appendix A. 

In the BCC structure Si represents overlap between the neighbors at distance (•\/3/2)i?, whereas for FCC it is 
{\/2/2)R, where R is the lattice parameters. We use these expressions to determine the self-adjusted wave functions 

To proceed with the numerical calculations a proper choice of the single-particle attractive potential in Eq. ^ has 
to be made which expressed the periodic kation potential energy. For that purpose, we have selected either 20 or 6 
potential wells surrounding given central ion, as shown in Fig. [21 

For such defined single-particle basis, the parameters Sa and t acquire the form for a linear chain 



Ea = P'To - ip-fTi + 2Y{To + T2), 
-t = t,,,+i - /32Ti - 2/37(To + T2) + 7'(3Ti + T3), 



(29a) 
(29b) 



where 



TfeEE(*,|iJi|*j) for jeSkii), 



(30) 



are the hopping integrals in the starting (Slater or Gaussian) basis. The corresponding quantities for the remaining 
lattices are provided in Appendix B. 

The corresponding expressions for U and Kij are more involved and will not be reproduced here. As the interaction 
parameters involve 6- fold integrals of the product of four wave functions, we have decided to use the Gaussian 
representation of ^i(r) of adjustable size rather than the Slater orbitals p^ iiii2,i3 ^ Namely, we shall use STO-nG 
basis. The coefficients of Gaussians are determined from their best adjustment to the Is Slater wave function. In 
Fig. [3] we exhibit the hydrogen Is state representation in representations ST0-3G, ST0-5G and ST0-7G, for which 
the ground state energy is respectively equal to —0.99169 i??/, — 0.99912 i?y and — 0.99987 i?y. Important is that the 
tail part playing the most important role in the narrow-band limit is reproduced quite well in the all cases. Also, one 
has to remember that the interaction parameters have a relatively large interatomic part. 
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FIG. 3: Comparison of Is Slater wave function (solid line) with its STO-3G, STO-5G and STO-7G representations. Inset: the 
difference in values near |r| — 0. 

IV. LIEB-WU SOLUTION COMBINED WITH AB INITIO ADJUSTMENT OF WANNIER 

FUNCTIONS 

A. Ground state properties 

The Hubbard-chain case (D = 1) with the exact Lieb-Wu (LW) diagonalization provides the canonical example of 
the application of EDABI method to the infinite systems. Here we treat also the LW solution as the only rigorous 
test of our approach. 

We start with the Lieb-Wu expressions including additionally e'^;^^ , as we will study the system evolution as a 
function of interatomic distance (not only as a function of model parameters, as they are varying with R). Explicitly, 
the ground-state energy expression is thenii 



Eg 



.eff 



At 



Jq{uj)Ji{uj) 



1 



pw(7/(4t) 



dw, 



(31) 



where J„(x) is the n-th order Bessel function. For the sake of comparison, the corresponding expression for the 
ground-state energy in the GWF approximation takes the form 



Eg 

N 



l"-4t / dkcosik)nk{g) + Ud{g), 



(32) 



where g is the Gutzwiller variational parameter, rik is the momentum distribution (seei^*^) and d expresses the 
probability of double occupancy. Also, the corresponding analytical expression in the Gutzwiller ansatz (GA) is 



^G_e/.f 4t/ nU/t 

N " n \ 32 



(33) 



These formulas represent a final step in solving model in LW, GWF and GA schemes, respectively. Here, these 
expressions represent a starting point for optimization of the single-particle basis {wi{r)}. We proceed as follows. 
First, we fix the interatomic distance R and for that distance construct the Wannier basis using the form of the periodic 
potential shown in Fig. [5] and the form P^ of the tight binding wave function. Then, we utilize the expression (PT]) 
to calculate Eg for given e^-^^ , t, U, which depend on the value of R. To be able to calculate the parameters U, Kij 
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FIG. 4: Ground state energy for the correlated chain vs. interatomic distance R (in units of Bohr radius ao) for the three 
methods of solving the extended Hubbard model, discussed in main text. Inset: R dependence of U/t ratio. 

(as well as t and e'^^^ in a single approach), we utilize the Gaussian basis ST0-3G and ST0-7G. In effect, the exact 
Lieb-Wu expression pip turns into a function, which is minimized with respect to a. This means that we change 
a ^ a + 6a and calculated Eq again and repeat the whole procedure until we reach the EQ{a) minimum for given 
R. This minimal energy represents the true physical ground-state energy and the corresponding value of a = amin 
as the inverse size of atomic wave-function in the correlated state. 

In Fig. m we display the ground state energy as a function of interatomic distance for LW, GWF and GA methods 
of approach. 

As one can expect, the exact (LW) solution provides the lowest energy for all values of R. However, the energy is 
higher than that of separated hydrogen atoms. So, the (extended) Hubbard chain with Is states is not stable in the 
vacuum. Such a quantum wire can be thus created only on an appropriate substrate. The atomic limit is practically 
reached for R > 5.5 ao, which corresponds to U/t > 50 (see the inset in Fig. [4]). For GA method the function Eg{R) 
ends at i? « 3.3 oq « 1.7 A for an optimized Wannier function with a « 1.021ao, /3 « 1.067, and 7 « 0.157. The 
important feature of the solution is that the U/t dependence on R is nonlinear, and thus those two representations of 
the results (either as a function of U/t or vs. R) are not equivalent. It is also important to note that in the inset in 
Fig. [4] we show U/t for all methods discussed; the curves practically coincide. Similar coincidence takes place also in 
the case of the parameters e'^;^-^, t, U, K = Ki^i+i, and related quantities, as shown in TableHl The same conclusion 
holds even when ST0-3G basis is utilized. One can see also that the atomic value oi U — (5/4) Ry is reached for 
i? ~ 6 flQ. However, the three methods provide quite different values of double occupancy probability d, as displayed 
in Fig. [3 



B. Discussion of the results for the Hubbard chain 



We first discuss the question to what extent the assumptions influence the numerical results. The first of them 
is the choice of the number of potential wells (cf. Fig. [5]). In Fig. [S]the ground state energy for LW solution using 
ST0-7G basis has been shown as a function of 1/fc, where k is the fc-th coordination sphere included in the procedure 
described in Sec. IIIIBI For fc > 6 we have already very reliable results. Second, the nearest-neighbor overlap S, 
justifying the tight-binding approximation, has been shown in Fig. [T] We see that for R/ao > 2 the overlap S is 
< 1/2. The difference between the ST0-3G and ST0-7G is small in that respect. These factors determine reliability 
of our results. Finally, in Fig. [5] we have plotted the nearest neighboring optimized Wannier functions when taking 
ST0-7G Gaussian basis. The characteristic cusps appearing at the lattice position have been smeared out to some 
extent. One should note that the functions contract by up 30% for a small lattice parameter. This is the reason why 
the inclusion of the electronic correlations extends the regime of applicability of the tight-binding approximation. One 



TABLE I: Microscopic parameters for the Hubbard chain as a function of lattice spacing for the three methods utilized. The 
ST0-7G basis was used. For details sec main text. 
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"The value for {R/ao)c ~ 3.288 is the threshold for the spurious metal-insulator transition within GA method. 
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FIG. 6: Ground state energy of the infinite Hubbard chain for the lattice parameter R — 2ao (ao - Bohr radius) obtained starting 
from the Lieb-Wu solution, plotted as a function of inverse k, where k denotes the number of coordination spheres taken to 
represent the potential V'(r). The limit 1/k — > represents the exact value for this lattice in the tight-binding approximation 
with Is orbitals. 
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FIG. 7: Overlap integral Si vs. interatomic distance, calculated for optimized the Wannier functions within the three schemes 
(LW, GWF, GA) and for the two selections of bases respectively. 



very important feature in Fig.[8]should be emphasized. Namely, the wave functions are almost the same independently 
of the method chosen to calculate them. This similarity explains why the microscopic parameters listed in Table I 
are close in value. Hence, it is the correlation energy which differentiates between the methods. 

One should also mention that the GWF solution requires determination of the statistical distribution function 
nk{g) {nk{g) = Uk] = JT-fcj, for paramagnetic, half- filled band case), which is contained in expression (|32p . This 
question has been analyzed in detail elsewherei^. The distribution function differs from that obtained exactly for 
nanoscopic systems-. Namely, in the latter case the statistical distribution n^ evolves from a Fermi-Dirac-momentum 
to continuous momentum distribution, without a characteristic cusp at Fermi momentum discussed by Metzner and 
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the Wannier functions are as follows: a = 1.194, /3 = 1.256, 7 = 0.341 (LW); a = 1.199,/? = 1.252, 7 = 0.338 (GWF); and 
a = 1.204, (5 = 1.248, 7 - 0.335 (GA). 



VoUhardO^. 

Strictly speaking, the expressions ([?T|) - ([55|) for the ground state energy are functionals of the wave functions {wi(r)}. 
So, to obtain the rcnormalized wave functions one should solve the Eulcr- Schrodinger equation which takes the form 



5Er 



Sw*{r) 



V- 



SEr 



,5V<(r) 



Xwi{r), 



(34) 



where A is the renormalizcd single-particle energy. This equation is very difhcult to solve when e.g. (|18[) is taken for 
Eg- This is a nonlinear, nonlocal differential- integral equation. To simplify the situation, we have chosen the trial 
wave functions with the variational parameter being the size of the selected Slater or Gaussian starting atomic wave 
function. 

In the next two sections we apply the same type of approach to higher-dimension lattices. Only GA is considered 
as the other two approaches are possible in the D = 1 case only. 



TWO DIMENSIONS 



For the lattices of the dimension D > 1 there is no exact analytical solution of Hubbard model, except for the 
solution in £) = 00 limit^^. Therefore, to illustrate our method for extended systems for 00 > D > 1 we use the GA 
solutio n^^i^^'^" . Within this scheme, the ground state energy can be calculated explicit and has the form 



Eg 

N 



.'^ff 



1 



U 



(35) 



where Uc — S\e\ is the critical value, for which we have a mean-field metal-insulator transition and e is the average 
bare band energy. The e is determined here from the expression, which in the paramagnetic case amounts to 



£ = iV-i(vI/o| 5] ty-a+ a,, 1*0) = - J2 ^W<0 



|k|<fef 



with 



e(k) 



Y^ 



t;,e 






ik(Ri-Rj) 



(36) 



(37) 



12 



TABLE 11: Average bare band energies (per site) for the lattices analyzed in this paper for the half-filled band. 



linear chain (CH) 
square lattice (SQ) 
triangular lattice (TR) 
simple cubic (SC) 
body center cubic (BCC) 
face center cubic (FCC) 



e = -(4/7r)t 
£ = -1.62 t 
e = -1.98 t 
e = -2.01 1 
e = -2.16 t 
£ = -2.64 i 



TABLE III: Selected microscopic parameters and quantities vs. lattice spacing for SQ and TR lattices within ST0-3G basis, 
using the GA expression for the ground state energy. The critical distance Rc/ao for the mean field metal-insulator transition 
is also specified. 





SQ 


TR 


SQ 


TR 


SQ 


TR 


SQ 


TR 


SQ 


TR 


R/ao 




(Ry) 


-t 


;%) 


U (Ry) 


K { 


Ry) 


Eg (Ry) 


2.0 


-0.135 


-0.034 


-0.494 


-0.357 


1.938 


2.087 


0.939 


0.946 


-0.524 


-0.315 


2.5 


-0.606 


-0.549 


-0.285 


-0.219 


1.701 


1.784 


0.769 


0.770 


-0.741 


-0.651 


3.0 


-0.808 


-0.777 


-0.181 


-0.144 


1.527 


1.579 


0.650 


0.649 


-0.844 


-0.804 


3.5 


-0.902 


-0.884 


-0.121 


-0.100 


1.418 


1.446 


0.563 


0.562 


-0.904 


-0.886 


iR/ao)c 


-0.906" 


-0.936' 


-0.090" 


-0.072' 


1.398" 


1.416' 


0.542" 


0.538' 


-0.919" 


-0.906' 


R/ao 


(3 


7 


a 


S 


d 


2.0 


1.417 


1.322 


0.295 


0.186 


1.233 


1.384 


0.466 


0.395 


0.174 


0.158 


2.5 


1.221 


1.207 


0.205 


0.146 


1.174 


1.241 


0.358 


0.323 


0.135 


0.121 


3.0 


1.129 


1.135 


0.153 


0.117 


1.107 


1.139 


0.281 


0.264 


0.087 


0.077 


3.5 


1.071 


1.083 


0.111 


0.091 


1.067 


1.078 


0.212 


0.207 


0.023 


0.022 


iR/ao)c 


1.057" 


1.068* 


0.099" 


0.082' 


1.063" 


1.067' 


0.192" 


0.187' 


0.000" 


o.ooo' 



"The value {R/ao)c ^ 3.648 for SQ. 
''The value (i?/ao)c ^ 3.668 for TR. 



being the bare band energy of particle with quasimonientuni k. In thcrmodynamical limit, s can be calculated through 
the integral 



e = 2/V* / e(k)d''k = 2 



dep{e)e, 



(38) 



Vf. 



e(k) 



where V* is the volume of the primitive unit cell in k space, and Vp represents the volume encompassed by the bare 
Fermi surface Sp. The second integral represents the integration over the density of states per atom spin (DOS) p{e) 
and with ep being the Fermi energy 

The numerical method of determining DOS was that of Buchheit and Loly^^. In the original work, the authors 
divided band into 100 intervals and the number of generated k was 3 • 10^. In this work, we divide the band into 500 
intervals with the 10^ generated wave vectors, which reproduce well the known analytic results for the chain. The 
average bare band energy per site obtained with the help of thus generated DOS is displayed for selected lattices in 
Table m 

In this Section we consider the square (SQ) and the triangular (TR) lattices. The single-particle attractive Coulomb 
potential was constructed for j g S'fc(i)U5fc(J + l) composed of fc = 13 spheres of coordination for SQ and k = 11 for TR 
cases respectively. In both situations such choice corresponds to up to five lattice parameters and involves 92 atomic 
sites for SQ and 102 atomic sites for TR cases. The comparison of the R dependence of selected microscopic parameters 
is provided in Table Hill Critical lattice parameter for metal-insulator transition (MIT) is Re ~ 3.648 oq ~ 1.930 A 
for SQ and Re ~ 3.985 ag ~ 2.109 A for TR. We remind the reader, that the corresponding quantity for the Hubbard 
chain in GA is Re ~ 3.318 oq ~ 1.756 A. The differences are caused by the different n.n number z = 2,4 and 6, 
respectively. 

The ground state energy Eq/N vs. R for both lattices is shown in Fig. [9l We see similarity for TR and SQ lattices, 
even though the U/t ratio similar for CH and SQ (cf. inset). The differences are included by different dependence of 
the double occupancy, as shown in Fig. 1101 The probability d is analytical function of U/t, namely 
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FIG. 9: Ground state energy vs. R for linear chain (CH) square lattice (SQ), and triangular lattice (TR). Inset: U jt ratio as 
a function of lattice parameter. The solid circles mark the MIT instability within GA solution. 
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The interesting quantity in this case is also the size a = l/a of the single atomic wave function composing the 
Wannier function. This quantity is plotted as a function of lattice parameter in Fig. [TlJ In all cases, correlations 
diminish the size of the starting atomic-like wave function. The size renormalization is 25-40% reduction for small 
lattice parameter. This result is counterintuitive, as one would expect, that the Coulomb repulsion increase the atomic 
size. This is not the case and, parenthetically, also the reason why the He atom, not H, is the smallest atom in the 
Universe^^. 
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FIG. 11: Inverse orbital size a vs. R of the atomic-like (Gaussian) wave function for the lattices specified. Inset: Same as a 
function of U/t. Note that the intersite interation is calculated in the Hartree-Fock approximation in the metallic state. 



VI. THREE DIMENSIONS 



In this Section we discuss the three cubic lattices. As in the case of 2D lattices, we select the sites generating the 
single-particle potential surrounding the central atomic site and one of its neighbors by atomic potentials distant by 
up to five lattice parameter. This case corresponds to the set of sites j S Sk{i) U Sk{i + 1) with fc = 22 (596 lattice 
sites) for SC structure, fc = 35 (1026 sites) for BCC structure, and fc = 47 (2344 sites) for FCC lattice. 

Selected parameters of Hubbard Hamiltonian are shown in Table HVl Note, that now R is the distance between the 
nearest neighbors, not the lattice parameter. 

In Fig. [13] we present the ground state energy vs. R As before, the Hubbard cubic lattice of Is are not stable on 
absolute energy scale. The double occupancy probability displayed in Fig. [TH The MIT transition at T = (from 
the paramagnetic metal to paramagnetic insulator) is continuous. For the sake of completeness, in Fig. [15] we have 
also shown the inverse size of atomic wave function vs. nearest-neighbor distance R. 

In Figs.[T2]we draw the nearest neighboring Wannier functions for SC structure along the directions [100] (left) and 
[111] (right) locating one of then at the origin. We see that overlap (and hopping) integral are strongly dependent on 
the direction, as one would expect. A similar effect is expected for the case of square lattice for D = 2. 



VII. DISCUSSION AND OUTLOOK 



In this paper we have shown on the example of a simple workable scheme - the Gutzwiller ansatz for the extended 
Hubbard model - how the calculation of single-particle wave function in the correlated state can be incorporated 
within the theory of correlated fermions in dimmensions D = 1,2,3. What is even more important, the dependence 
of the physical properties (e.g. the ground state energy) is calculated as a function of interatomic distance and this 
dependence is not related linearly to that as a function of the parameter U/t. In other words, only the explicit 
dependence on R provides the proper scaling of the quantities. Additionally, the Wannier orbitals calculated from 
the variational principle for the energy of the correlated-electron state are renormalized in the sense that they are 
more tightly bound to their parent ions than their bare correspondants, thus extending the applicability of the tight- 
binding approach for narrow-band systems. These are the most important qualitative conclusions coining from our 
method of approach. In essence, as for the correlated system the interaction is comparable (or even stronger) than 
the single-particle part, the method reverts the usual approach, i.e. we determine the single-particle wave function 
only after diagonalization of the many-particle parametrized Hamiltonian. Such a method is fully implemented here, 
albeit for the case of a narrow s-band. 
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TABLE IV: Selected microscopic parameters for cubic lattices vs. R. 
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FIG. 12: The n.n. renormalized Wannier functions along [100] direction for SC structure (left); right: the functions along [111] 
direction. 



The importance of the Hubbard model in its present form for the real 3d systems such as oxides, manganites or 
intermetallic compounds, is limited. This is because the orbital degenerancy of the 3d states is not included explicitly 
and the interorbital Hunds's rule and direct Coulomb interactions do not appear here. Therefore, the factors such 
as the Hund's rule are neglected. It would be important to extend there present method by including crystal- field 
effects due to anions and interatomic interorbital interactions. The first steps in this direction have been undertaken 
alreadyS^. The same remark concerns the inclusion of the 3d-2p hybrydization in the metallic oxides and other related 
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FIG. 13: Ground state energy vs. the nearest-neighbor distance R for the cubic lattices. Inset: I] jt as a function of R for those 
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FIG. 14: Double occupancy probability d = {ni^riii) as a function of the distance R between nearest neighbors for the cubic 
lattices. Inset: Double occupancy probability vs. U/t. 



compounds. On the example of our model calculations we see that while the single-particle wave function does not 
depend strongly on the method of treating the many-particle Hamiltonian in the Fock space, (cf. Sec. IVB), the 
two-particle correlation functions (e.g. d{R)), are strongly method-dependent. This is the reason why many-body 
renormalizations of the quasiparticle characteristics are not computable in a reliable manner within any of the present 
band-structure calculations. So, extensions of this method to the degenerate Hubbard model solution combined with 
wave-function adjustment in the correlated state, is highly desirable. 
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FIG. 15: Inverse size a vs. R of the atomic-like (Gaussian) wave function for the lattices specified. Inset: Same as a function 
of U/t. 
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APPENDIX A: THE EXPRESSIONS FOR COEFFICIENTS OF A AND B 

Below we list the explicit expressions for the coefficients A and B defined in Scc. lIIIBJ for the linear chain (CH) 

A=l + S2, B = 3 + S3, (Ala) 

the square lattice (SQ) 

A=1 + 2S2 + S3, B = 9 + ASi + Se, (Alb) 

the triangular lattice (TR) 



the simple cubic (SC) 



the body-centered cubic (BCC) 



A = I + 2S1 + 2S2 + S3, 

B = l + 15S'i + 6S'2 + 6S'3 + 6S'4 + S'5, 



A^l + 4S2 + Si, B^m + SSa + US^ + Ss, 



A ^ I + 3S2 + 3S3 + S5, 

B = 27S1 + 27S4 + 9S7 + S10 



(Ale) 
(Aid) 

(Ale) 
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and the face-centered cubic (FCC) 



A = l+4S'i + 252 + 4S'3 + 5'4, 

B = A + AbSi + 1252 + 365.3 + 1254 + 1855 + 456 + 1257 + S'g- (Alf) 

These lattices coefficients are used in the subsequent analysis. 

APPENDIX B: THE EXPRESSIONS FOR COEFFICIENTS OF £„ AND t FOR THE ALL LATTICES 

UNDER CONSIDERATION 

For the single-particle basis defined in subsection IIIIBl the parameters £a and t acquire the form for a linear chain 

£a = /?2ro-4/37ri + 272(ro + r2), (Bla) 

-t = i,,,+i = l3^Ti - 2l3-f{To + T2) + j^STi + T3). (Bib) 

For the square lattice (SQ) the single parameters are expressed through their atomic (Gaussian) basis correspondents 
as follows 

Ea = P^To~SfijTi+A-r^{To + 2T2+T3), (B2a) 

-t = p^Ti-2p-f{To + 2T2 + T3)+-/^{9Ti + 6Ti + Te). (B2b) 

For the triangular lattice they are 

Ea = l3To - 12f3jTi + 67^(^0 + 2Ti + 2r2 + T3), (B3a) 

-t = p^Ti-2p-f{To + 2Ti + 2T2+T3)+j\2To + l5Ti + 6T2 + 6T3 + 6T^ + T5). (B3b) 

For the cubic lattices the single particle parameters are expressed through their atomic (Gaussian) basis correspon- 
dents as follows: 
for SC 

Ea = p^To-12(3jTi + 6^^{To + AT2+Ti), (B4a) 

-t = p^Ti-2p-f{To + 4T2+Ti)+j^{15Ti + 8T3 + 12T5+Ts), (B4b) 



for BCC 



and for FCC 



Sa = l3''To-l6p-fTi+8-f\To + 3T2 + 3T3 + n), (B5a) 

-t = f3''Ti-2Pj{To + 3T2 + 'ST3 + T5)+j\27Ti + 27T^ + 9TT + Tio), (B5b) 



Sa = /32T0 -24/37ri + 127^(^0 + 4T1 + 2T2+ 473 +T4), (B6a) 

-t = /32T1 - 2/37(To + 4ri + 2T2 + 4T3 + 74) + 

7^ {4To + 45Ti + I2T2 + 3673 + 1274 + ISTg + 476 + I2T7 + Tg). (B6b) 
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^ For review see: V.I. Anisimov, J. Zaanen, and O.K. Anderson, Phys. Rev. B 44, 993 (1991); V.I. Anisimov, F. Aryasetiawan, 

and A.I. Lichtenstein, J. Phys.: Condens. Matter, 9, 767 (1997). 
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